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We present a comparison between finite differences schemes and a pseudospectral method applied to the 
numerical integration of stochastic partial differential equations that model surface growth. We have studied, in 
1+1 dimensions, the Kardar, Parisi and Zhang model (KPZ) and the Lai, Das Sarma and Villain model (LDV). 
The pseudospectral method appears to be the most stable for a given time step for both models. This means that 
the time up to which we can follow the temporal evolution of a given system is larger for the pseudospectral 

■ method. Moreover, for the KPZ model, a pseudospectral scheme gives results closer to the predictions of the 
O ' continuum model than those obtained through finite difference methods. On the other hand, some numerical 
D , instabilities appearing with finite difference methods for the LDV model are absent when a pseudospectral 

integration is performed. These numerical instabilities give rise to an approximate multiscaling observed in 
the numerical simulations. With the pseudospectral approach no multiscaling is seen in agreement with the 
^ continuum model. 

> 

C/3 , PACS numbers; 81. 15.Aa,05.40.-a,64.60.Ht,05.70.Ln 

■ I. INTRODUCTION 

o 

Q ' Kinetic surface roughening of surfaces growing in nonequilibrium conditions has been intensively studied for the last two 
I— ' decades 1110. Htl- Theoretical approaches make use of both discrete atomistic simulations and stochastic continuum equations 
for the evolution of the coarse-grained surface height /i(x, t). There is overwhelming experimental evidence that surfaces under 
, general nonequilibrium growth conditions can develop scale-invariant correlations in space and time, which supports the hope of 
' a unified theoretical framework to understand kinetic roughening phenomena from first principles. The aim is at identifying the 
\^ various dynamical universalities of growth associated with different sets of symmetries and/or conservation laws. It is believed 
that only these basic elements largely determine the universality class and the value of the coiTesponding critical exponents. In 

■ theoretical studies attention is therefore focused on symmetries and only the most relevant terms (in the renormalization group 
sense) are expected to be required to describe a particular class of growth. 

Universality classes of growth are generically represented by stochastic partial differential equations, 

i> : 

O' dth^g{vh) + T^{^,t), (1) 

. y—i, where /i(x, t) is the height of the interface at substrate position x and time t. The external noise ?7(x, t) represents the influx of 
atoms on the surface. The function Q{Vh) defines a particular model and incorporates the relevant symmetries and conservation 
5h laws. In particular, invariance under translation along the growth and substrate directions as well as invariance in the election 

. . . 1 of the time origin rule out an explicit dependence of Q on h, x and t. Very often the presence of nonlinearities in Q require 
the use of perturbative renormalization techniques to obtain analytical approximations for the critical exponents, which can then 
be compared with Monte Carlo simulations of atomistic models and experiments. A perturbative renormalization approach 
invariably provides the critical exponents as a series expansion on the parameter e ~ dc — d, where the critical dimension dc 
can be very high when compared with the dimensions of physical interest (usually d = 1 or 2). Only in a few lucky cases some 
extra symmetries produce cancellation of higher order loop diagrams that results in a scaling relation between exponents to be 
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exact to all orders in the perturbative expansion. More often than not, the case is that we only have approximations to the critical 
exponents valid up to a certain order in e and a great deal of elaborated algebraic effort is required to improve our approximation 
up to the next order This often makes direct numerical integration of Eq. ([TJ an extremely useful and necessary tool as the most 
reliable source of precise values for the critical exponents. 

Numerical schemes to integrate continuum surface growth equations like Eq. ([T]l in 1 + 1 and 2 + 1 dimensions tend to be 
unsophisticated. In most cases a straightforward finite-differences (FD) method on a lattice does an excellent job and provides 
highly precise values for the critical exponents, at least in dimensions of experimental relevance. In this approach (see details 
below) one basically approximates the continuous height field, h{x,t), by its values on the lattice sites, hi{t), and derivatives 
by differences between neighboring sites. More clever choices of the discretization rule have been shown to be useful to obtain 
better agreement with exact properties of the continuum solutions |4], which could not be obtained by using a conventional 
discretization, like the nominal values of the continuum equation parameters. 

However, the use of FD schemes sometimes poses some important problems fHI^Sl- In particular Dasgupta et. al. have 
shown !9t] by means of numerical simulations that discretized versions of commonly studied nonlinear growth equations 
exhibit a instability in the sense that single pillars (grooves) become unstable when their height (depth) exceeds a critical value. 
In some cases these instabilities are not present in the corresponding continuum equations, indicating that the behavior of the 
discretized versions is indeed different from their continuum counterparts. It is important to remark that this pillar/groove 
instability is actually generic to the FD discretizations of a large class of nonline ar g rowth equations, including the Kardar- 
Parisi-Zhang (KPZ) [10] or the Lai-Das Sarma- Villain (LDV) equations HI [H [Tj . This is a puzzling result because the 
corresponding continuum equations are not unstable. 

In many situations, like for instance in KPZ, the existence of this instability is of little significance for practical purposes and 
one can actually carry out a correct numerical integration by using FD schemes. The reason is that one is mostly interested 
in the growth from a flat (or almost) surface initial condition and common relaxation mechanisms do not favor the formation 
of large pillars or grooves. In these cases the instability is only realized if the initial condition is prepared in such a state that 
there is a pillar/groove of size above the threshold on a otherwise flat surface, which is highly artificial and usually uninteresting 
for practical purposes. However, as already pointed out in Ref. [9], there is a large class of systems for which the instability 
of any FD scheme is inevitable. Specifically, discrete versions of models exhibiting anomalous kinetic roug hening HSH 
will certainly show this kind of instability at sufficiently long times. The reason being that anomalous scaling is 
associated with a nontrivial dynamics of the average surface gradient (local slope), so that ((V/i)^)^/^ ^ i", with k > lfl4l[l7ll . 
Therefore, systems exhibiting anomalous roughening will dynamically generate large local height differences, no matter how flat 
the initial condition is. As a consequence, provided that a simulation is run long enough, the surface will produce pillars/grooves 
above the critical value for the instability to appear, like for instance is the case for LDV. 

One of the models we study in this paper is the LDV equation ifTll [T2I [Tsll 

/i(x, t) = -KV^h + AV2(V/i)2 + ^(x, t), (2) 
where the noise is Gaussian distributed and delta correlated, 

(e(x, t)e(x', t')) = 2D5{^ - ^)6{t - t'). (3) 

This model constitutes a minimal model for the long wavelength behavior of surface growth under ideal molecular-beam epitaxy 
conditions. The LDV model is interesting in many respects and has been the focus of a lot of attention in the literature iflU \Va 

[H mils El HI- 

Numerical simulations of discrete versions of (|2]i in 1 + 1 dimensions have reported 11911 a finite, albeit small, anomalous 
exponent k k. 0.08, possibly indicating a logarithmic dependence. A theoretical prediction fir] based on Flory-type arguments 
predicted n « 1/11 (see however |23]). Therefore, from the discussion above, one would expect a discrete version of (|2]i to 
become unstable. This problem was studied by Dasgupta et. al. and they showed ijsllStl that FD algorithms were actually unstable 
at long times. They also estimated the critical height step to be around hc{\) « A/\ with A « 20 for dU with K — D = \, 
which clearly shows that the instability will appear the sooner the larger the nonlinear coefficient is. Those authors also found 
that the addition of higher-order nonlinearities in the FD version of the model controls the numerical instabilities and renders 
a stable surface, but with intermittent fluctuations and multiscaling properties of surface correlations. It has been claimed that 
higher-order nonlinearities of the form A2,i (Vft.)^", with n > 1, may play an important role in LDV universality class because 
they are infinitely many marginally relevant nonlinear terms 1*7^, [Tu [Til I20I1 . 

These results can be compared with FD integration schemes for the KPZ equation 

(9f/i(x,i) = i^V2/i + A(V/i)2 +e(x,t). (4) 

It has been shown f^,^ that discrete versions of Eq. (|4|l were stable, unless isolated grooves of large enough size are included 
in the initial state. The reason being that KPZ exhibits conventional (no anomalous) scaling and local slopes are thus rapidly 
converging towards a constant. Under general conditions the constant is much smaller than the critical slope he for the KPZ 
discretization to become unstable and so, large slopes are not spontaneously generated by the dynamics. 
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In reference ['29'] a study of the ID and 2D KPZ equation using finite-differences and pseudospectral integration methods 
is presented. Authors claim that a pseudospectral method gives results closer to the continuum limit than finite-differences 
methods. They show how a pseudospectral method reproduces the exact value of the global width of the steady state interface 
within error bars whereas a finite-differences method with conventional discretization of the nonlinear term entails significant 
differences in the amplitude value. They also use pseudospectral computations to reproduce the most reliable values of the 
dynamical exponents obtained through discrete growth models. 

In this paper we discuss the validity of FD integration algorithms in the presence of anomalous roughening. We compare 
the accuracy, stability and overall performance of FD methods versus pseudospectral (PS) schemes applied to the paradigmatic 
examples of the KPZ and LDV equations. We claim that the instability previously found in FD discretizations is spurious and 
non physical, therefore, FD should be generally avoided in numerical simulations of continuum growth models with anomalous 
scaling. We argue that the main reason for the adequacy of PS to attack growth problems with anomalous scaling is that spatial 
derivatives are more accurate than in FD methods, where one implicitly assumes that the step height is small. Our conclusions 
are based on numerical analysis of KPZ and LDV equations in 1 + 1 dimensions by means of FD and PS integration schemes. 
Our results when comparing both techniques are conclusive: (i) PS methods are stable against isolated pillars/grooves, while 
FD are not, ( ii) under the same conditions PS schemes take much longer than FD to get to a numerical overflow, and ( Hi) PS 
schemes give well behaved correlation functions with no trace of multiscaling. Finally, we will discuss the implications of our 
results for the appearance of multiscahng in discrete models proved to be in the same universality class as LDV equations ll24ll . 



II. NUMERICAL INTEGRATION SCHEMES 



In order to perform a numerical integration of Eqs. (|4]l and (|2]l the parameters can easily be rescaled to have only one inde- 
pendent control parameter- namely, the coupling constant. As it is customary, one can work with dimensionless variables h, x 
and t so that all parameters but one are set to unity, so we have 

dfh{x,t) = \7^h + g{Vhf +r]{x,t), (5) 

for the KPZ equation, where the dimensionless coupling constant is g ~ \\J1D jv^. We can also write the LDV equation in 
dimensionless form 

dth(y., t) = -V^/i + g^'^iyhf + ^^(x, i), (6) 

where the coupling constant is g — \^2D / and 77 (x, t) is a Gaussian noise with mean zero, unit variance and correlations 
(?7(x, i)ry(x', t')) ^ 5(x - ^')6{t - t'). 

Let us now summarize the idea behind FD and PS integration schemes and introduce some useful definitions. Equations (ID 
and (|6]l can be cast in the form 

a* /i(x, t)=C [h] (x, i) + $ [h] (x, t) + 77(x, t) , (7) 
where C[h] is a linear functional of h and is another functional containing the nonlinear terms. 



A. Finite-differences metliods 



We consider a rf-dimensional lattice with periodic boundary conditions with uniform spacing Ax in each direction. The 
positions of the nodes in the lattice are given by 

Xj - Aa;(ji, j2, . . . , Jd), ^ sC iV, - 1, 1 ^ ^ < d. (8) 

where iV, is the lattice size in the i-th direction. Using a one step Euler's method to compute the temporal derivative, the 
evolution of a system governed by Eq. (|7]i is given by: 



/i(xj, t + M) - /i(xj, t) + At(/:[/^](xj, t) + a>H(xj, t)) + ^ rKxj, i). (9) 

where At is the time step and the stochastic variables 77(xj, t) have zero mean and correlations (?7(xj, i)77(xj, t')) = Sj yS{t~t'). 
We took the i] variables as Gaussian random numbers (other distributions can be used as long as they satisfy the central hmit 
theorem). 
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In finite difference methods, derivatives are computed by truncating the Taylor series of the field up to certain order Let us 
introduce the finite difference operators and which are, respectively, the forward and backward difference operators along 
the j direction: 

A5/i(x, t) = /i(x + ejAx) - /i(x, t), (10) 
A]hiyi,t)=h{yi,t)-h{ii-ejAx). (11) 
In terms of these operators, the linear parts of Eqs. (|5]l and (|6]l are, up to second order of approximation, given by: 

d 

£KPz(xj,t) = (v2/i)(xj,i) = (Ax)-2^A^A^/^(xj,i), 

i=l 

/:LDv(xj,<) = -(V4/i)(xj,i) = -V\V'h) ^ {Axr' f2 At'A^A5A5/i(xj,i). 

The explicit expressions in 1 + 1 dimensions are 

£kpzN = iAx)-^h,+i - 2h, + 
CLDv[h] = -iAx)-^{hi+2 - 4/ii+i + 6h, - 4/i,_i + 

where Xi = iAx, i = 0, . . . , N ~ 1 are the positions of the nodes in the lattice and hi ~ h{xi, t). Regarding the nonlinear terms 
we consider for the gradient square the usual symmetric discretization: 

(VM'(x,i) = \{Ax)-^Y.[{A\ + A\)h{^,t)f 

1=1 

that in 1 + 1 dimensions becomes 

{Vhf{x,,t) = ]^{Ax)-\h,+^ ^ h,^if . (12) 

In the case of the KPZ equation, other discretizations of the nonlinear term have been proposed ll4l l25ll . We mention the Lam 
and Shin discretization (LS) 10] 

(V/i)2(xj,t) = \{Ax)-^Y.{[{A\ + A\)h{^;,t)f - (AfMxj,t))(A^Mxj,t))}, 

i=i 

so that in 1+1 dimensions we have 

{Vhf{x,,t) = i(Aa;)-2[(/i,+i - + - K){h, - + {h, - h.^^f]. (13) 

LS discretization has two interesting features in 1+1 dimensions: (i) the effective parameter g agrees with its nominal value, 
and ( ii) the probability distribution of the discretized version in the steady state can be computed exactly and it turns out to be 
the probability distribution of the continuum equation for all values of g. It has been argued that this discretization allows 
to recover some results predicted by the continuum model while discrepancies when using the conventional discretization (fT2l i 
have been observed |4]. 

In the following we use a lattice spacing Ax — 1. As is customary in this kind of simulations, hydrodynamic limit is achieved 
by increasing the number of lattice sites N . In numerical integrations of continuous growth models one avoids to perform the 
Ax — > limit with fixed L, which would lead the system towards the linear critical point, since the coupling constant of the 
discretized equation is 5 ^ as Ax ^ 0. A fixed lattice spacing Ax in the limit L ^ 00 is always preferred as it best drives 
the system towards the nontrivial critical point. 

B. Pseudospectral method 

To compare with FD methods we have considered a numerical scheme consisting of a spectral method in space together with a 
Euler's method in time. We assume that the field ft,(x, t) satisfies periodic boundary conditions in the multidimensional interval 
[0, LY and we represent it as a truncated Fourier series 

ker„ 
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The set Tn over which the sum is taken is given by Fat = {(fci, fc2, . . . , kd) / ~ N/2 ^ fc^ ^ N/2 — 1, 1 i ^ d}, and the 
hk{tys are the Fourier coefficients of h, defined as 



^ J\O.L]'i 



The noise term 77 is also replaced by its expansion rj^ in Fourier modes. When N 00 the usual Fourier series is recovered. 
On the other hand, when h and 77 are replaced by and rj^ respectively in Eq. Q, the residual 

i?Ar(x, t) = dthN - C[hN] - ^[Hn] - Vn 

will be not null in general. By requiring Rn to be orthogonal to the functions {e^"^'^, k G Fat}, we obtain a set of ODEs for the 
Fourier coefficients of h. This procedure is actually equivalent to project the equation onto a subspace of orthogonal polynomials 
of degree ^ N/2. Then, by imposing the orthogonality condition 



/ dxi?Ar(x,i)e""»'' = 0, keFjv 



we obtain 

dt 



= cjk/ik(<) + *k(t)+?7k(t), kGFjv. (14) 



The quantity is the linear dispersion relation, which is obtained through the Fourier transform of the linear part of the 
equation; it is Wk = <^ for Eq. ^ and Wk = — for Eq. (|6]l. The <E>k(i)'s are the Fourier coefficients of the nonlinear terms 
and are given by the following convolution sums: 



-9 X! ■ ^2 Ki K2 (KPZ), 

ki+k2=k 

9^ X! qi • ^2 /iki ^k2 (LDV). 

ki+k2=k 



Regarding the Fourier coefficients of the noise, it is easy to verify that ?7k(i) are complex Gaussian variables with zero mean and 
correlations 

{fi^.{t)fn.'{t')) = ^5^^^^,5{t~t'). 

The Fourier coefficients /ik are in general difficult to compute. In addition, even the simplest nonlinearities make it computa- 
tionally expensive the task of going from real space to Fourier space and viceversa. For these reasons, we consider a discretized 
space with N nodes in each direction 

= ^0'i'i2,-.-,jd), o^j, i-^i-^d, 

and we use the discrete Fourier transform to integrate ( fT4l i. The discrete Fourier coefficients depend only on the values of the 
field at the nodes xj and are given by (direct discrete Fourier transform): 



j 

where hj{t) — /i(xj, t). We have the inversion formula (inverse discrete Fourier transform): 



kerj. 



Then, we replace the continuum Fourier coefficients in ( fT4b by the discrete ones, so that 



^Wk/jkW + ^kW + 'ykW- (15) 
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Note that Eq. (fTsT i is now written in terms of the discrete Fourier coefficients /ik- To integrate ( fTSl ) we perform the following 
change of variables based on the solution of the linear equation: 

/ik(i) = e^'^'z^it) + i?k(t) 

where 

Jo 

The Zk's satisfy the equations 



dt 



^kWe"""*- (16) 



The set of ODEs dTSI l can be solved by using one of the several algorithms available for stochastic differential equations (Euler, 
Runge-Kutta, predictor-corrector methods, etc.). Considering a one step Euler's method to integrate ( fTSI l and going back to the 
original variable /ik, we are finally left with 

/ik(i + AO =e"''^*[/ik(t) + At4k(i)]+fk(t). (17) 

Equation ( fTTj i can be reinterpreted as an Euler scheme with time-step (the factor multiplying the nonlinear term) e^^^^/St, so 
our algorithm provides a smaller time-step for the smallest length-scales (so it is intrinsically multiscale). This represents a 
significant improvement with respect to the pseudospectral method used in li29tl which is just Eq. ( fTSl l integrated with a one step 
Euler's method. 

Assuming that Wk = ^-k, the variables fk {t) can be obtained as 



'"-^^ = V 2c.k (A.)^ ^"-(^^ 

where Ax = L/N and v\^{t) are the discrete Fourier transform of a set of Gaussian random numbers of zero mean and unit 
variance. Note that, as expected, when At ^ we recover the last term on the right side of Eq. (|9|l. 

The computation of the nonlinear terms for the KPZ and LDV equations in Fourier space involves the Fourier transform of 
the product of two functions (actually, the square of V/i). In general, calling these two functions a and p, we need to calculate 
the convolution sum 



*k= '7kiPk2- (18) 



ki+k2=k 

ki,k2eriv 

In one dimension, this convolution sum implies 0{N'^) operations, which is computationally more expensive than a finite 
difference method, for which only 0{N) operations are needed. To speed up the computation, we used a pseudospectral 
transform method to compute the Fourier transform of the nonlinear term. Starting from (Tk and pk, the inverse transformation is 
used to obtain a and p in real space. Then a and p are multiplied to obtain $ in real space. Finally, the direct Fourier transform 
is applied to obtain the $k- In terms of the discrete Fourier operator JF, this pseudospectral calculation can be written as follows: 

This procedure allows to evaluate the convolution sum using 0{N log N) operations in one dimension. It is important to note that 
the Fourier coefficients $k computed in a pseudospectral manner differ from those obtained from a true spectral computation. 
The difference is the so-called aliasing error. For example, in one dimension, the coefficients $fe computed pseudospectrally 
turn out to be 

ki+k2 = k ki+k2=k±N 

The first term on the right hand side is just the convolution sum ( fTSl l whereas the second term is the aliasing error The aliasing 
error has been proved to be asymptotically of the same order of the error made in truncating the Fourier series. There are several 
recipes to remove the aliasing. We used a well-known truncation technique usually referred to as the 3/2-rule li27ll . 
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FIG. 1 : Average over space and realizations of the nonlinear term for tlie ID KPZ equation as a function of time. The inset shows the average 
height of the interface as a function of time. Here L — 128, and the averages are talcen over 100 realizations. The kind of numerical method 
and the value of the nonlinear coupling parameter for each curve are shown in the legend. 



III. COMPARISON OF THE METHODS IN 1+1 DIMENSIONS 



A. Preliminaries 

In order to compare the results provided by the FD and PS numerical methods applied to models ^ and we must first 
notice that, for any given model and the same value of the nonlinear coupling parameter g, the intensity of the nonlinear effects 
depends on the numerical scheme used to integrate the equation. This fact, which has already been pointed out in [28], leads 
to the conclusion that different algorithms cannot in principle be compared directly. In Fig.[T] the average (both over space and 
reaUzations) of the nonlinear term for the ID KPZ equation is shown in several cases. We see that, for the same value of the 
coupling parameter, the PS method gives effectively a larger nonlinear term. In other words, for the same value of g, the FD 
method underestimates the intensity of the nonlinear term with respect to the PS method. A comparison of the two numerical 
methods can be only made if the nonlinear effects are of the same order for both of them on average. 

For the KPZ equation the nonlinear effects can be monitored by measuring the mean velocity of the interface, which is given 

by 

v=^ dx{{Vhf). 
L Jo 

In the inset of Fig. [T]the average height of the interface as a function of time for the ID KPZ equation is shown. The slope of 
this curve is just the velocity of the interface. For the same value of g, the interface obtained with the PS method moves faster 
As said before, this indicates that the nonlinear effects are stronger in the PS method than in the FD method. It is easy to find 
values of g such that the interface in both cases moves approximately at the same velocity, which means that nonlinear effects 
are of similar magnitude. Then, if we denote by wfd(5) and Vps{g) the mean interface velocity for the FD and PS methods, 
respectively, the value of the coupling parameter g such that wps(ff) = wfd(5) is given by 

^FD(g) 

vps{9) 

The ratio WFD(3)/wps(g) depends smoothly on both g and the system size as it can be seen in Fig. |2] The ratio of velocities 
of the interfaces slightly decreases with g and increases with the system size. For example, for a system size of L = 128 and 
g = 2.5, numerically we have found that vpu{g)/vps{g) ~ 0.46, which means that the nonlinear effects in the PS method are 
approximately twice as much stronger than those of the FD method. In Fig. [1] we can see how the nonlinear terms for both 
integration methods become similar when g is decreased from 2.5 to a value of 2.5 x 0.46 = 1.15 for the PS method. 

In the case of the LDV equation, we can proceed in a similar manner. Let us denote by ipM{g',t) — g{\W^{\/h)'^\) the absolute 
value of the nonlinear term of the LDV equation averaged over space and realizations for the numerical method M. Then, for a 
given value of g used with the FD method, we can estimate a g for the PS method leading to a nonlinear term of the same order. 
This is 

\V'ps(5;0/t 
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FIG. 2: Ratio ufd/wps as explained in the text as a function of the nonhnear couphng parameter g for several system sizes. 




FIG. 3: Average over space and realizations of the absolute value of the nonlinear term for the ID LDV equation as a function of time. Here 
L — 128, and the averages are taken over 100 realizations. The kind of numerical method and the value of the nonlinear coupling parameter 
for each curve are shown in the legend. 



In the previous expression the angular brackets denote an average over the time interval used in the simulation. As for the KPZ 
equation, the ratio g/g computed according to Eq. ( |20b depends slightly on both the system size and g. For example, for a 
system size of L = 128 and a value of gpD — 1-25 used in the FD algorithm, we find that a value of gps ~ 0.42 gives rise to 
a nonlinear term of the same order for the PS method (see Fig. [3]). As occurs for the KPZ equation, the nonlinear effects are 
stronger for the PS method. 

We checked that the global dynamical exponents obtained with the FD and PS methods are the same using values of g 
according to ( fT9l ) and (l20b . The global interface width scales according to the Family- Vicsek ansatz jsoll : 

WiL,t) {{h{x,t)-hY) ^ <"/^/(L/<i/^), 
where the scaling function / behaves as 




It" ifu<l, 
const if M 3> 1. 



The parameter a is the roughness exponent, z is the dynamic exponent, and the ratio (3 = a/z is the time exponent. In 1+1 
dimensions the critical exponents can be computed exactly 1 10] and their values are a = 1/2 and z = 3/2, so that /3 = 1/3. 
Using the FD with g — 2.5, we found the exponents a ~ 0.49, j3 ~ 0.32, z — a/ (3 ~ 1.52 for the FD method and with g = 1.2 
we obtain the same values of the exponents with the PS method within error bars. For the LDV equation the global exponents 
are known for arbitrary dimension. In 1+1 dimensions they are a = 1, /3 = 1/3, and z — 2>. Taking a value of g — 1.25 we 
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FIG. 4: Probability of instability for the ID KPZ equation as a function of time. The initial condition is a flat interface. Here L = 100 and 
the probabilities are computed using 2000 realizations. Curves for two values of the time step are shown. We show results for the FD and PS 
numerical methods and some different values of the nonlinear coupling parameter g which are shown in the legend. 



found the exact value of the exponents with two significant digits by integrating the equation with the FD method, with system 
sizes ranging from L = 16 to L = 256 and averaging the interfaces over 100 runs. On the other hand, the PS method with 
g = 0.42 provides the same exponents within error bars. 



B. Stability of the algorithms 

We tested the stability of the algorithms by measuring the probability P{t) that the system exhibits a numerical overflow when 
starting from a flat interface. This is measured from a large number of independent runs as the frequency probability of getting 
a computer overflow at time t. This numerical instability takes place when the height of the interface tends to grow indefinitely. 
The probability of instability is a decreasing function of the time step used in the simulations. 

In Fig. m the probability of instability as a function of time for the ID KPZ equation is shown for several cases. We show 
curves for two time steps At = 10~^ and At = 10~^. The system size is L = 100 and the probabilities are computed over 2000 
samples. The values of g were chosen in such a way that the nonlinear effects for the two methods were of the same order. For 
the KPZ equation this is achieved when gps ~ 0.48(7fd. In all cases we found the probabilities for the PS method to be smaller 
than those of the FD method for a given time step. For example, for a time step of Ai = 10^'^, we can see in the bottom graphic 
of Fig.|4]that the PS method is stable (that is, the probability of instability is equal to zero) in the time interval [0, lOO] for values 
of g — 3.6, 4.8, and 6.0, whereas the FD becomes unstable at very short times. 

In Fig.|5]we show the probability of instability as a function of time for the ID LDV equation. In this case we took ^ps — 
0.34(7fd to match the nonlinear effects for both methods. In much the same way as for the KPZ equation we see that for a given 
time step the PS method is the most stable. This is also observed for other time steps ranging from 10"^ to 10"^ We then 
conclude that the PS method is the most stable when the intensity of the nonlinear terms are of equivalent magnitude. This 
means that, under the same conditions, the PS method allows to follow the temporal evolution of the system up to larger times. 



C. KPZ equation 



There are some exact results of the continuum KPZ model that we can used to test the numerical methods. First, the steady 
state probability distribution of the heights is known exactly |[lll2l- In terms of the slopes, m^x) — dxh{x), it is known that 



P{m) ~ exp 
This expression can be written approximately as: 

/ „ 

V{m) ~ exp [ " A, 



dx m[xY 



exp (— Aa; m 



1=0 
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FIG. 5: Probability of instability for the ID LDV equation as a function of time. The initial condition is a flat interface. Here L = 100 and 
the probabilities are computed using 2000 realizations. Curves for two values of the time step are shown. We show results for the FD and PS 
numerical methods and some different values of the nonlinear coupling parameter g which are shown in the legend. 



We have used the central limit theorem to identify X]r=o ^"i ^^'^^ Here m represents the slope of the field in the steady state 
at any point of the lattice. The normalized expression of the probability is: 

r{m) = 7r-i/2g-m^^ 

With this probability distribution, we can find the several moments of the slope to: {m?) — 1 /2, {ra^) — 3/4, (to^) — 15/8. 
For the discretized model and when the LS discretization is used, the steady state distribution probability is found to be 101 



P[hi\ ^ Gxp 



N 



(21) 



It is worth mentioning here a caveat concerning ( 1211 1. One can see that, in order to reproduce numerically the slope distribution 
with the FD method, slopes must be computed with the forward (or backward) operator (Eq. (fTOli). so that mi = h[ = hi+i — hi. 
We have checked that if the symmetric rule to compute the derivatives {rrii — (/i^+i — /ii_i)/2) is used the width of the slope 
probability distribution is far from unity, which is the exact value. When forward or backward derivatives are used, however, the 
correct value is recovered. Remarkably, the PS method provides the proper result in a natural way. 

The global interface width in the steady state is also known exactly 1I21I1 . 

W{L) = ^J^^L^'\ t^^, (22) 

and it is independent of the nonlinear coupling parameter g. In reference it is shown that a FD method with conventional 
discretization for the nonlinear term, Eq. ( fT2] i. provides steady state interfaces whose global width is of the form ( l22b but 
with a prefactor of V-/"^ significantly smaller than the predicted value 24"^^^. It has been argued |4] that with the improved 
discretization ( TOt the correct value for the prefactor is recovered. A plot of 4){L) — y/2A/LW{L) versus L^^ was presented in 
references i23.l29ll . showing that 4>{L) is unity within error bars for both the PS method and the FD method with the discretization 
( fT3] l. although the dispersion of the data is larger for the FD method. We have also carried out a similar study comparing FD and 
PS methods. We checked that the curve W{L) vs. L can be fitted to a function of the form B L^/"^, where B = 0.182 ± 0.002 
for the FD method with the usual discretization ( fT2] i. As expected, this value is clearly smaller than the nominal value Bq = 
24-1/2 ^ 0.204. This observation is in agreement with that of reference [4]. On the other hand, for both PS and FD method 
with LS improved discretization (fT3] ) we obtain the same value (indistinguishable up to the third digit) B = 0.196 ± 0.003, a 
value very close to Bq indeed. Therefore, with the PS method we obtain in a natural way the result predicted by the continuum 
model for the steady state global interface width. For the FD method, on the contrary, we must use an ad hoc discretization of 
the nonlinear terms to achieve the same results. 
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FIG. 6: Ratios aq/ai, q = 2,3, 4, 5, as a function of time, of tlie moments of the nearest neighbor height difference for the ID LDV equation 
integrated with the FD [Fig. (a)] and PS [Fig. (b)] numerical methods. Here L = 100, (?fd ~ 2, gps — 1.5, and the averages are taken over 
300 realizations. 



D. LDV equation 

We have investigated the influence of the numerical method on the multiscaling behaviour of the LD V eq uation. The multi- 
scahng can be detected by looking at the moments of the nearest neighbour height difference. We define (]M 

'J,{t) = {\K+l{t)-K{t)\'y^': 

where the average is taken over each site of the system and over the different realizations of the noise. In systems that exhibit 
anomalous scaling, as the case of the LDV equation for intermediate times (see note ll23ll ). the moments aq{t) are expected to 
grow as follows 

where L is the system size and z is the dynamical exponent. When the moments scale in a different way, that is, when the ag's 
depend on g, the system is said to show multiscaling. As done in lUllst], we monitorize the multiscaling by looking at the ratios 
crq{t)/aiit),q > 2. On the other hand, the multiscaling can also be studied by measuring the height difference correlation 
function |laS[3- In Fig. |6^ we show the ratios aq{t) / ai{t) with q = 2,3,4,5 for the ID LDV equation integrated with 
the FD scheme. Parameter values are L = 1000, g — 2, and the averages are taken over 100 realizations. This is in fact a 
reproduction of Fig. 12 of reference As can be seen in the picture, the greater the q, the faster the growth of aq/tJi with 
time. This behaviour implies the existence of multiscaling. For times greater than 100, the evolution of the system cannot be 
followed due to the presence of numerical instabilities. The authors of pi] claimed that this behavior is related to an instability of 
the discretized LDV equation against the growth of isolated pillars, which are just height profiles such that the field h is positive 
at a certain point while being zero otherwise. For a given a value of the parameter g, there exists a critical value he of the pillar 
height beyond which the pillar grows with a certain probability. It is found that this critical height goes as he ~ g^^- On the 
other hand, the effect of the magnitude of the time step on this instability seems to be very small. For this reason, it is further 
argued in [9] that this instability is not a numerical artifact due to the use of a not small enough integration time step. As we 
show in the following our results disagree with this interpretation. 

Interestingly, we find that the numerical integration of the LDV equation with the PS method has a significant impact on the 
observed scaling properties. The instability discussed above, which is present with any FD method, is no longer present, at least 
in the wide range of couplings we studied. Specifically, a piUar like initial condition of the form 



otherwise. 



with s$ j TV - 1, TV odd, /iq > 



(23) 
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never grows with the PS method. In the absence of noise, the temporal derivative of the field ( l23l l is always negative. A 
straightforward calculation leads to 

d^h = -V^h + gV\Vhf = -^AN^ - l)(3iV2 - 7) « < 0, L = TV » 1. 

Therefore, a pillar or the form ( |23]) always tends to shrink in the deterministic case. Other structures Uke double pillars, however, 
might grow in time when their size exceeds a certain value. 

In Fig. |6}3 we show the ratios aq (t) /cti (t), q ~ 2,3,4,5, for the ID LDV equation integrated with the PS scheme. In this case 
we use a value of g — 1.5 for which the nonlinearities are considerably stronger than for the FD method with g — 2. As it can 
be observed, the curves do not grow in time, which means that there is no multiscaling. For other values of g and other system 
sizes similar results were obtained. It is worth mentioning, however, that the instability does show up when the aliasing is not 
removed. In this case an isolated pillar may grow when its height is larger than a critical value, in much the same way as with 
the FD method. So, it is strongly recommended to remove aliasing effects when applying a PS method to correctly describe the 
continuum physics. Note that the aliasing does not actually affect global properties of the system such as multiscaling, the value 
of the critical exponents, etc. Indeed, we have also carried out a study of the PS method without removing the aliasing because 
in this case a similar instability of the FD method is present. We have not observed multiscaling in this case either We then 
conclude that the multiscaling does not seem to be related to this instability for the PS method. 

This analysis leads to the main conclusion of our paper that the existence of multiscaling may depend on the numerical scheme 
used to integrate the growth model. We argue that the instability (and associated multiscaling behavior) is intrinsic to the nu- 
merical integration scheme rather than to the discretization itself, in contrast with the conclusions of Ref. jsi |^. Our results 
clearly show that the instability previously found in FD discretizations has to be seen as spurious and inherent to the discretiza- 
tion scheme used. PS integration methods do not show any trace of either instability (when aliasing is properly removed) or 
multiscaling, representing much more accurately the dynamics and statistics of the continuum problem. We conclude that a PS 
method should be preferred for surface growth equations with anomalous scaling, because next-neighbor height differences in 
this case can grow very large. 

IV. CONCLUSIONS 

We have shown that the choice of the numerical method used to integrate certain stochastic models of surface growth may be 
of paramount importance in the study of some physical properties of the system. We have compared a standard finite difference 
method with a pseudospectral numerical scheme in the integration of the KPZ and LDV growth models in 1h-1 dimensions. As 
the FD method underestimates the nonlinear effects with respect to the PS method, the nonlinear coupling parameter were tuned 
up so that the nonlinear terms were of the same order for both numerical methods on average. The global critical exponents, 
obtained from the global interface width are the same for the two numerical methods. 

With regard to the KPZ equation there are some exact results available derived from the continuum model. The expression 
for the global width of the interface is known in the saturation regime. With the FD method and a standard discretization for the 
nonlinear term, the amplitude of the width of the numerical interfaces is smaller than that of the continuum. With the spectral 
approach, on the contrary, numerical results are very close to the predicted value. Nevertheless, it is possible with the finite 
differences scheme to get close to the continuum model prediction for the width of steady state interfaces, but at the expense of 
using more sophisticated discretizations. 

We have tested the stability of the algorithms by measuring for different time steps the probability of the system to undergo 
a floating point instability evolving from a flat interface. This instability is related to a numerical overflow in the surface height 
data. The PS method proved to be the most stable in all the cases we have studied for both models. In the same way, with the PS 
method it is possible to follow the temporal evolution of the system for longer times than with the FD method. 

The LDV equation exhibits anomalous scaling at intermediate times, so that according to 1 17] (but see also li23ll ) the average 
slope of the field is expected to grow in time. Any FD method leads to a numerical instability against the growth of an isolated 
pillar appears. This instability has been claimed |8, 9] to be the reason why approximate multiscaling is observed in the numerical 
simulations, although multiscaling is not present in the continuum equation. More importantly, this multiscaling has been 
interpreted as a real physical effect, which could explain the multiscaling behavior of surface fluctuations observed in atomistic 
models believed to belong to the LDV universality class. However, our results show that this interpretation is misleading. We 
have shown that surface multiscaling is not observed with the PS method, regardless of the temporal evolution of isolated pillars. 
Therefore, surface multiscaling does not seem to be related to this instability for the PS method nor represent any intrinsic physics 
of the LDV equation as such. In this respect, due to the fact that discrete models can be mapped into continuum equations, in 
particular the LDV equation f2?], our results indicate that multiscaling behavior observed in such systems could be an artifact 
of the discretization of the dynamics and, consequently, the are not intrinsic to the physical system they are trying to modehze. 
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